#check means and variances of all variables
#varTable(F2fit)
#model summary
#summary(F2fit, fit.measures=TRUE)
#compare the two models
#cbind(m1=lavInspect(F1fit, 'fit.measures'), m2=lavInspect(F2fit, 'fit.measures'))
#anova(F1fit, F2fit)
#Perceived Discrimination and Peer Victimization Scales
Discrimdata <- select(fulldata, "safe1", "safe2", "safe8", "safe10", "safe13", "safe15", "safe16", "safe20", "safe21",
"bully4", "bully5", "bully6", "bully7")
#look over summary statistics for items
summary(Discrimdata)
str(Discrimdata)
#recheck summary statistics
summary(Discrimdata)
str(Discrimdata)
apply(Discrimdata, 2, table)
#fit one factor model
F1.model <- 'discbully =~ safe1 +safe2 +safe8 + safe10 + safe13 + safe15 + safe16 + safe20 + safe21 + bully4 + bully5 + bully 6 + bully7'
F1fit <- cfa(F1.model, data=Discrimdata, missing = "fiml") #using fiml to handle missing data
#check means and variances of all variables
#varTable(F1fit)
#model summary
summary(F1fit, fit.measures=TRUE)
#fit two factor model
disc.model <- 'safe =~ safe1 +safe2 +safe8 + safe10 + safe13 + safe15 + safe16 + safe20 + safe21
bully =~ bully4 + bully5 + bully 6 + bully7'
discfit <- cfa(disc.model, data=Discrimdata, missing = "fiml") #using fiml to handle missing data
#check means and variances of all variables
varTable(discfit)
#model summary
summary(discfit, fit.measures=TRUE, standardize = TRUE)
#compare the two models
#cbind(m1=lavInspect(F1fit, 'fit.measures'), m2=lavInspect(F2fit, 'fit.measures'))
anova(F1fit, discfit)
#InternalizingSympotoms Factor Model
intnames <- paste(rep("ysr_prob_", times = 32), c(14, 29, 30, 31, 32, 33, 35, 45, 50, 52, 71, 91, 112, 5, 42, 65, 69, 75, 102, 103, 111, 47, 51, 54, "56a", "56b", "56c", "56d", "56e", "56f", "56g", "56h"), sep = "")
intdata <- select(fulldata, intnames)
#look over summary statistics for items
summary(intdata)
str(intdata)
apply(intdata, 2, table)
#fit one factor model
intform <- paste(intnames, collapse = "+")
int.model <- paste("internalizing =~ ", intform)
intfit <- cfa(int.model, data=intdata, missing = "fiml") #using fiml to handle missing data
#check means and variances of all variables
varTable(intfit)
#model summary
summary(intfit, fit.measures=TRUE, standardized = TRUE)
View(fulldata)
fulldata$safescale <- fulldata$safe1 +fulldata$safe2 +fulldata$safe8 + fulldata$safe10 + fulldata$safe13 + fulldata$safe15 + fulldata$safe16 + fulldata$safe20 + fulldata$safe21
fulldata$bullyscale <- fulldata$bully4 + fulldata$bully5 + fulldata$bully6 + fulldata$bully7
fulldata$familismoscale <- fulldata$FS1 + fulldata$FS2 + fulldata$FS3 + fulldata$FS4 + fulldata$FS5 + fulldata$FS6 + fulldata$FS7
fulldata$internalizingscale <- fulldata$ysr_prob_14+fulldata$ysr_prob_29+fulldata$ysr_prob_30+fulldata$ysr_prob_31+fulldata$ysr_prob_32+fulldata$ysr_prob_33+fulldata$ysr_prob_35+fulldata$ysr_prob_45+fulldata$ysr_prob_50+fulldata$ysr_prob_52+fulldata$ysr_prob_71+fulldata$ysr_prob_91+fulldata$ysr_prob_112+fulldata$ysr_prob_5+fulldata$ysr_prob_42+fulldata$ysr_prob_65+fulldata$ysr_prob_69+fulldata$ysr_prob_75+fulldata$ysr_prob_102+fulldata$ysr_prob_103+fulldata$ysr_prob_111+fulldata$ysr_prob_47+fulldata$ysr_prob_51+fulldata$ysr_prob_54+fulldata$ysr_prob_56a+fulldata$ysr_prob_56b+fulldata$ysr_prob_56c+fulldata$ysr_prob_56d+fulldata$ysr_prob_56e+fulldata$ysr_prob_56f+fulldata$ysr_prob_56g+fulldata$ysr_prob_56h
View(fulldata)
fulldata$safescale
table(fulldata$Dem1a)
table(fulldata$Dem2)
fulldata$anxdepscale <- fulldata$ysr_prob_14+fulldata$ysr_prob_29+fulldata$ysr_prob_30+fulldata$ysr_prob_31+fulldata$ysr_prob_32+fulldata$ysr_prob_33+fulldata$ysr_prob_35+fulldata$ysr_prob_45+fulldata$ysr_prob_50+fulldata$ysr_prob_52+fulldata$ysr_prob_71+fulldata$ysr_prob_91+fulldata$ysr_prob_112
fulldata$withdepscale <- fulldata$ysr_prob_5+fulldata$ysr_prob_42+fulldata$ysr_prob_65+fulldata$ysr_prob_69+fulldata$ysr_prob_75+fulldata$ysr_prob_102+fulldata$ysr_prob_103+fulldata$ysr_prob_111
fulldata$somaticscale <- fulldata$ysr_prob_47+fulldata$ysr_prob_51+fulldata$ysr_prob_54+fulldata$ysr_prob_56a+fulldata$ysr_prob_56b+fulldata$ysr_prob_56c+fulldata$ysr_prob_56d+fulldata$ysr_prob_56e+fulldata$ysr_prob_56f+fulldata$ysr_prob_56g+fulldata$ysr_prob_56h
fulldata$anxdepscale
Dem1a
table(fulldata$Dem1a)
dat2keep <- cbind(fulldata$Dem1a,fulldata$Dem2,fulldata$safescale,fulldata$bullyscale,fulldata$familismoscale,fulldata$internalizingscale)
dat2keep
dat2keep <- cbind(seq(1,nrow(fulldata),by = 1),fulldata$Dem1a,fulldata$Dem2,fulldata$safescale,fulldata$bullyscale,fulldata$familismoscale,fulldata$internalizingscale)
dat2keep
dat2keep[is.na] <- 999
dat2keep[is.na,] <- 999
dat2keep[is.na(dat2keep)] <- 999
dat2keep
write.table(dat2keep, "Ramos.dat", col.names = F, row.names = F)
table(fulldata$Dem2)
fulldata$Dem2 <- fulldata$Dem2 - 1
View(dat2keep)
fdir::set()
fdir::set()
fdir::set()
dat <- read.table("imp1.dat")
names <- c("id","age","gender","discrim","bully","familism","internal","female")
names(dat) <- names
View(dat)
fdir::set()
dat <- read.table("imp1.dat")
names <- c("id","age","gender","discrim","bully","familism","internal","female")
names(dat) <- names
summary(dat[c("discrim","bully","familism","internal")])
#Factor models for Predictors of Internalizing Symptoms in Rural Latinx Youth: A Serial Mediation Model of Sociocultural Factors
#install.packages("dplyr")
#install.packages("lavaan")
#libraries
library(dplyr) #used for select function
library(lavaan) #used for cfa function
fdir::set()
#read in data
olddata <- "FINAL 3 Levels  Dummy Code Serial Mediation with update safe syntax 2.19.20 DD.csv"
newdata <- "FINAL Serial Mediation_4.3.20.csv"
fulldata <- read.csv(newdata, header = TRUE)
fulldata$Dem2 <- as.factor(fulldata$Dem2)
levels(fulldata$Dem2) <- c("Male", "Female")
fulldata <- filter(fulldata, between(Dem1a, 14, 17))
#555 is missing for bully items
fulldata$bully4[fulldata$bully4==555] <- NA
fulldata$bully5[fulldata$bully5==555] <- NA
fulldata$bully6[fulldata$bully6==555] <- NA
fulldata$bully7[fulldata$bully7==555] <- NA
#Non-integers are missing for safe items
fulldata$safe1[fulldata$safe1%%1 != 0] <- NA
fulldata$safe2[fulldata$safe2%%1 != 0] <- NA
fulldata$safe8[fulldata$safe8%%1 != 0] <- NA
fulldata$safe10[fulldata$safe10%%1 != 0] <- NA
fulldata$safe13[fulldata$safe13%%1 != 0] <- NA
fulldata$safe15[fulldata$safe15%%1 != 0] <- NA
fulldata$safe16[fulldata$safe16%%1 != 0] <- NA
fulldata$safe20[fulldata$safe20%%1 != 0] <- NA
fulldata$safe21[fulldata$safe21%%1 != 0] <- NA
fulldata$bully5[fulldata$bully5%%1 != 0] <- NA
fulldata$bully6[fulldata$bully6%%1 != 0] <- NA
#Familismo and family conflict scales
#Note: Family conflict was removed, and only familismo is used in the final model
#names(fulldata)
familydata <- select(fulldata,"FS1","FS2","FS3","FS4","FS5","FS6","FS7")
#look over summary statistics for items
summary(familydata)
str(familydata)
#some added tables to look at frequency of observed cases
apply(familydata, 2, table)
#fixing values for scales (ie., FA item scores can only go from 1-5, fri scores can only be 0 or 1)
#itemdata$FS2[itemdata$FS2==0] <- 1
#itemdata$FS4[itemdata$FS4==0] <- 1
#itemdata$FS5[itemdata$FS5==0] <- 1
#itemdata$FS6[itemdata$FS6==0] <- 1
#itemdata$FS7[itemdata$FS7==0] <- 1
#itemdata$fri3[itemdata$fri3==30] <- NA
#recheck summary statistics
summary(familydata)
str(familydata)
#fit one factor model
Family.model <- 'familismo =~ FS1 + FS2 + FS3 + FS4 + FS5 + FS6 + FS7'
#using fiml to handle missing data
Familyfit <- cfa(Family.model, data=familydata)
summary(Familyfit, fit.measures = TRUE)
#check means and variances of all variables
varTable(Familyfit)
#model summary
summary(Familyfit, fit.measures=TRUE)
#fit two factor model
#F2.model <- 'familismo =~ FS1 + FS2 + FS3 + FS4 + FS5 + FS6 + FS7
#             familyconflict =~ fri3 + fri6 + fri9 + fri12 + fri15 + fri18 + fri21 + fri24 + fri27'
#using fiml to handle missing data
#F2fit <- cfa(F2.model, data=itemdata,
#             ordered = c("fri3", "fri6", "fri9", "fri12", "fri15", "fri18", "fri21", "fri24", "fri27"))
#check means and variances of all variables
#varTable(F2fit)
#model summary
#summary(F2fit, fit.measures=TRUE)
#compare the two models
#cbind(m1=lavInspect(F1fit, 'fit.measures'), m2=lavInspect(F2fit, 'fit.measures'))
#anova(F1fit, F2fit)
#Perceived Discrimination and Peer Victimization Scales
Discrimdata <- select(fulldata, "safe1", "safe2", "safe8", "safe10", "safe13", "safe15", "safe16", "safe20", "safe21",
"bully4", "bully5", "bully6", "bully7")
#look over summary statistics for items
summary(Discrimdata)
str(Discrimdata)
#recheck summary statistics
summary(Discrimdata)
str(Discrimdata)
apply(Discrimdata, 2, table)
#fit one factor model
F1.model <- 'discbully =~ safe1 +safe2 +safe8 + safe10 + safe13 + safe15 + safe16 + safe20 + safe21 + bully4 + bully5 + bully 6 + bully7'
F1fit <- cfa(F1.model, data=Discrimdata, missing = "fiml") #using fiml to handle missing data
#check means and variances of all variables
#varTable(F1fit)
#model summary
summary(F1fit, fit.measures=TRUE)
#fit two factor model
disc.model <- 'safe =~ safe1 +safe2 +safe8 + safe10 + safe13 + safe15 + safe16 + safe20 + safe21
bully =~ bully4 + bully5 + bully 6 + bully7'
discfit <- cfa(disc.model, data=Discrimdata, missing = "fiml") #using fiml to handle missing data
#check means and variances of all variables
varTable(discfit)
#model summary
summary(discfit, fit.measures=TRUE, standardize = TRUE)
#compare the two models
#cbind(m1=lavInspect(F1fit, 'fit.measures'), m2=lavInspect(F2fit, 'fit.measures'))
anova(F1fit, discfit)
#InternalizingSympotoms Factor Model
intnames <- paste(rep("ysr_prob_", times = 32), c(14, 29, 30, 31, 32, 33, 35, 45, 50, 52, 71, 91, 112, 5, 42, 65, 69, 75, 102, 103, 111, 47, 51, 54, "56a", "56b", "56c", "56d", "56e", "56f", "56g", "56h"), sep = "")
intdata <- select(fulldata, intnames)
#look over summary statistics for items
summary(intdata)
str(intdata)
apply(intdata, 2, table)
#fit one factor model
intform <- paste(intnames, collapse = "+")
int.model <- paste("internalizing =~ ", intform)
intfit <- cfa(int.model, data=intdata, missing = "fiml") #using fiml to handle missing data
#check means and variances of all variables
varTable(intfit)
#model summary
summary(intfit, fit.measures=TRUE, standardized = TRUE)
fulldata$safescale <- fulldata$safe1 +fulldata$safe2 +fulldata$safe8 + fulldata$safe10 + fulldata$safe13 + fulldata$safe15 + fulldata$safe16 + fulldata$safe20 + fulldata$safe21
fulldata$bullyscale <- fulldata$bully4 + fulldata$bully5 + fulldata$bully6 + fulldata$bully7
fulldata$familismoscale <- fulldata$FS1 + fulldata$FS2 + fulldata$FS3 + fulldata$FS4 + fulldata$FS5 + fulldata$FS6 + fulldata$FS7
fulldata$internalizingscale <- fulldata$ysr_prob_14+fulldata$ysr_prob_29+fulldata$ysr_prob_30+fulldata$ysr_prob_31+fulldata$ysr_prob_32+fulldata$ysr_prob_33+fulldata$ysr_prob_35+fulldata$ysr_prob_45+fulldata$ysr_prob_50+fulldata$ysr_prob_52+fulldata$ysr_prob_71+fulldata$ysr_prob_91+fulldata$ysr_prob_112+fulldata$ysr_prob_5+fulldata$ysr_prob_42+fulldata$ysr_prob_65+fulldata$ysr_prob_69+fulldata$ysr_prob_75+fulldata$ysr_prob_102+fulldata$ysr_prob_103+fulldata$ysr_prob_111+fulldata$ysr_prob_47+fulldata$ysr_prob_51+fulldata$ysr_prob_54+fulldata$ysr_prob_56a+fulldata$ysr_prob_56b+fulldata$ysr_prob_56c+fulldata$ysr_prob_56d+fulldata$ysr_prob_56e+fulldata$ysr_prob_56f+fulldata$ysr_prob_56g+fulldata$ysr_prob_56h
fulldata$anxdepscale <- fulldata$ysr_prob_14+fulldata$ysr_prob_29+fulldata$ysr_prob_30+fulldata$ysr_prob_31+fulldata$ysr_prob_32+fulldata$ysr_prob_33+fulldata$ysr_prob_35+fulldata$ysr_prob_45+fulldata$ysr_prob_50+fulldata$ysr_prob_52+fulldata$ysr_prob_71+fulldata$ysr_prob_91+fulldata$ysr_prob_112
fulldata$withdepscale <- fulldata$ysr_prob_5+fulldata$ysr_prob_42+fulldata$ysr_prob_65+fulldata$ysr_prob_69+fulldata$ysr_prob_75+fulldata$ysr_prob_102+fulldata$ysr_prob_103+fulldata$ysr_prob_111
fulldata$somaticscale <- fulldata$ysr_prob_47+fulldata$ysr_prob_51+fulldata$ysr_prob_54+fulldata$ysr_prob_56a+fulldata$ysr_prob_56b+fulldata$ysr_prob_56c+fulldata$ysr_prob_56d+fulldata$ysr_prob_56e+fulldata$ysr_prob_56f+fulldata$ysr_prob_56g+fulldata$ysr_prob_56h
dat2keep <- cbind(seq(1,nrow(fulldata),by = 1),fulldata$Dem1a,fulldata$Dem2,fulldata$safescale,fulldata$bullyscale,fulldata$familismoscale,fulldata$internalizingscale)
dat2keep
names(dat2keep) <- names
summary(dat2keep)
for(v in 2:ncol(dat2keep)){
plot(density(dat2keep[,names[v]], na.rm = T), main = names[v])
}
summary(dat2keep)
v <- 2
plot(density(dat2keep[,names[v]], na.rm = T), main = names[v])
names <- c("id","age","gender","discrim","bully","familism","internal")
names[v]
plot(density(dat2keep[,names[v]], na.rm = T), main = names[v])
dat2keep[,names[v]]
names[v]
dat2keep[,v]
summary(dat2keep)
for(v in 2:ncol(dat2keep)){
plot(density(dat2keep[,v], na.rm = T), main = names[v])
}
summary(dat2keep)
names(dat2keep) <- names
summary(dat2keep)
names
dat2keep
dat2keep <- cbind(seq(1,nrow(fulldata),by = 1),fulldata$Dem1a,fulldata$Dem2,fulldata$safescale,fulldata$bullyscale,fulldata$familismoscale,fulldata$internalizingscale)
names <- c("id","age","gender","discrim","bully","familism","internal")
names(dat2keep) <- names
dat2keep
dat2keep <- cbind(seq(1,nrow(fulldata),by = 1),fulldata$Dem1a,fulldata$Dem2,fulldata$safescale,fulldata$bullyscale,fulldata$familismoscale,fulldata$internalizingscale)
names <- c("id","age","gender","discrim","bully","familism","internal")
dat2keep
colnames(dat2keep) <- names
dat2keep
summary(dat2keep)
fdir::set()
dat <- read.table("imp1.dat")
names <- c("id","age","gender","discrim","bully","familism","internal","female")
names(dat) <- names
summary(dat[c("discrim","bully","familism","internal")])
for(v in 2:ncol(dat)){
plot(density(dat[,v], na.rm = T), main = names[v])
}
vars.2.sim <- c("discrim","bully","familism","internal")
dat.2.sim <- dat[vars.2.sim]
for(v in 3:ncol(dat.2.sim)){
res <- rnorm(nrow(dat.2.sim), 0, 2)
dat.2.sim[, v] <- round(dat.2.sim[, v] + res, rnd[v])
}
vars.2.sim <- c("discrim","bully","familism","internal")
dat.2.sim <- dat[vars.2.sim]
for(v in 3:ncol(dat.2.sim)){
res <- rnorm(nrow(dat.2.sim), 0, 2)
dat.2.sim[, v] <- round(dat.2.sim[, v] + res, 0)
}
dat.2.sim
vars.2.sim <- c("discrim","bully","familism","internal")
dat.2.sim <- dat[vars.2.sim]
for(v in 1:ncol(dat.2.sim)){
res <- rnorm(nrow(dat.2.sim), 0, 2)
dat.2.sim[, v] <- round(dat.2.sim[, v] + res, 0)
}
dat.2.sim
vars.2.sim <- c("discrim","bully","familism","internal")
dat.2.sim <- dat[vars.2.sim]
for(v in 1:ncol(dat.2.sim)){
res <- rnorm(nrow(dat.2.sim), 0, 2)
dat.2.sim[, v] <- round(dat.2.sim[, v] + res, 0)
dat.2.sim[dat.2.sim[,v] <= 0, v] <- 0
}
dat.2.sim
cor(dat,dat.2.sim)
for(v in 2:ncol(dat.2.sim)){
plot(density(dat.2.sim[,v], na.rm = T), main = vars.2.sim[v])
}
dat.2.save <- cbind(dat$id, dat$age, dat$female, vars.2.sim)
vars.2.sim
dat.2.save <- cbind(dat$id, dat$age, dat$female, dat.2.sim)
dat.2.save
write.csv(dat.2.save, "ramos.csv")
# example x: moderated regression with fiml estimation and auxiliary variables
library(fdir)
library(lavaan)
library(mdmb)
##################################################
# read data
##################################################
# set working directory
set()
# read raw data from working directory
data <- read.table("behaviorachievement.dat", na.strings = "999")
names(data) <- c("id","male","hispanic","riskgrp","atrisk","behsymp1","lrnprob1",
"read1","read2","read3","read9","readgrp9","stanread7",
"math1","math2","math3","math9","mathgrp9","stanmath7")
model <- "stanread7 ~ 1; read2 ~ 1; hispanic ~ 1; read9 ~ 1; read1 ~ 1; lrnprob1 ~ 1; behsymp1 ~ 1;"
descriptives <- inspectSampleCov(model, data, meanstructure = T, missing = "fiml")
# center lower-order variables
data$read1.cgm <- data$read1 - descriptives$mean["read1"]
descriptives$mean["read1"]
data$lrnprob1.cgm <- data$lrnprob1 - descriptives$mean["lrnprob1"]
descriptives$mean["lrnprob1"]
descriptives <- inspectSampleCov(model, data, missing = "fiml")
descriptives$mean["read1"]
# summarize incomplete predictors to determine ranges for pseudo-imputations
summary(data[,c("stanread7","read2","hispanic","read9","read1.cgm","lrnprob1.cgm","behsymp1")])
# example x: moderated regression with fiml estimation and auxiliary variables
library(fdir)
library(lavaan)
library(mdmb)
##################################################
# read data
##################################################
# set working directory
set()
# read raw data from working directory
data <- read.table("behaviorachievement.dat", na.strings = "999")
names(data) <- c("id","male","hispanic","riskgrp","atrisk","behsymp1","lrnprob1",
"read1","read2","read3","read9","readgrp9","stanread7",
"math1","math2","math3","math9","mathgrp9","stanmath7")
##################################################
# estimate means and center variables
##################################################
# estimate sample statistics in lavaan
model <- "stanread7 ~ 1; read2 ~ 1; hispanic ~ 1; read9 ~ 1; read1 ~ 1; lrnprob1 ~ 1; behsymp1 ~ 1;"
descriptives <- inspectSampleCov(model, data, missing = "fiml")
# center lower-order variables
data$read1.cgm <- data$read1 - descriptives$mean["read1"]
data$lrnprob1.cgm <- data$lrnprob1 - descriptives$mean["lrnprob1"]
##################################################
# analysis
##################################################
# summarize incomplete predictors to determine ranges for pseudo-imputations
summary(data[,c("stanread7","read2","hispanic","read9","read1.cgm","lrnprob1.cgm","behsymp1")])
# set ranges (nodes) for pseudo-imputations
nodes.hispanic <- c(0,1)
nodes.read2 <- seq(10, 160, by = 5)
nodes.stanread7 <- seq(80, 420, by = 10)
nodes.read9 <- seq(30, 130, by = 2)
nodes.read1 <- seq(-55, 75, by = 2)
nodes.lrnprob1 <- seq(-30, 50, by = 2)
nodes.behsymp1 <- seq(10, 100, by = 2)
# model for behsymp1 predictor: f(behsymp1)
model.behsymp1 <- list( "model" = "linreg", "formula" = behsymp1 ~ 1, nodes = nodes.behsymp1)
# model for lrnprob1 predictor: f(lrnprob1|behsymp1)
model.lrnprob1 <- list( "model" = "linreg", "formula" = lrnprob1.cgm ~ behsymp1, nodes = nodes.lrnprob1)
# model for read1 predictor: f(read1|lrnprob1,behsymp1)
model.read1 <- list( "model" = "linreg", "formula" = read1.cgm ~ lrnprob1.cgm + behsymp1, nodes = nodes.read1)
# model for read9 outcome: f(read9|read1,lrnprob1,behsymp1)
model.read9 <- list( "model" = "linreg", "formula" = read9 ~ read1.cgm + lrnprob1.cgm + read1.cgm*lrnprob1.cgm + behsymp1, nodes = nodes.read9)
# model for hispanic auxiliary variable f(hispanic|read9,read1,lrnprob1,behsymp1)
model.hispanic <- list("model" = "logistic", "formula" = hispanic ~ read9 + read1.cgm + lrnprob1.cgm + behsymp1, nodes = nodes.hispanic)
# model for read2 auxiliary variable f(read2|hispanic,read9,read1,lrnprob1,behsymp1)
model.read2 <- list("model" = "linreg", "formula" = read2 ~ hispanic + read9 + read1.cgm + lrnprob1.cgm + behsymp1, nodes = nodes.read2)
# model for stanread7 auxiliary variable f(stanread7|read2,read9,read1,lrnprob1,behsymp1)
model.stanread7 <- list("model" = "linreg", "formula" = stanread7 ~ read2 + hispanic + read9 + read1.cgm + lrnprob1.cgm + behsymp1, nodes = nodes.stanread7)
# combine predictor models into a list
predictor.models <- list(read2 = model.read2, hispanic = model.hispanic, read9 = model.read9, behsymp1 = model.behsymp1, lrnprob1 = model.lrnprob1, read1 = model.read1)
# estimate factored regression model w mdmb
fit <- frm_em(dat = data, dep = model.stanread7, ind = predictor.models)
summary(fit)
# summarize incomplete predictors to determine ranges for pseudo-imputations
summary(data[,c("stanread7","read2","hispanic","read9","read1.cgm","lrnprob1.cgm","behsymp1")])
# set ranges (nodes) for pseudo-imputations
nodes.hispanic <- c(0,1)
nodes.read2 <- seq(10, 160, by = 5)
nodes.stanread7 <- seq(90, 410, by = 10)
nodes.read9 <- seq(30, 130, by = 2)
nodes.read1 <- seq(-55, 75, by = 2)
nodes.lrnprob1 <- seq(-30, 50, by = 2)
nodes.behsymp1 <- seq(10, 100, by = 2)
# model for behsymp1 predictor: f(behsymp1)
model.behsymp1 <- list( "model" = "linreg", "formula" = behsymp1 ~ 1, nodes = nodes.behsymp1)
# model for lrnprob1 predictor: f(lrnprob1|behsymp1)
model.lrnprob1 <- list( "model" = "linreg", "formula" = lrnprob1.cgm ~ behsymp1, nodes = nodes.lrnprob1)
# model for read1 predictor: f(read1|lrnprob1,behsymp1)
model.read1 <- list( "model" = "linreg", "formula" = read1.cgm ~ lrnprob1.cgm + behsymp1, nodes = nodes.read1)
# model for read9 outcome: f(read9|read1,lrnprob1,behsymp1)
model.read9 <- list( "model" = "linreg", "formula" = read9 ~ read1.cgm + lrnprob1.cgm + read1.cgm*lrnprob1.cgm + behsymp1, nodes = nodes.read9)
# model for hispanic auxiliary variable f(hispanic|read9,read1,lrnprob1,behsymp1)
model.hispanic <- list("model" = "logistic", "formula" = hispanic ~ read9 + read1.cgm + lrnprob1.cgm + behsymp1, nodes = nodes.hispanic)
# model for read2 auxiliary variable f(read2|hispanic,read9,read1,lrnprob1,behsymp1)
model.read2 <- list("model" = "linreg", "formula" = read2 ~ hispanic + read9 + read1.cgm + lrnprob1.cgm + behsymp1, nodes = nodes.read2)
# model for stanread7 auxiliary variable f(stanread7|read2,read9,read1,lrnprob1,behsymp1)
model.stanread7 <- list("model" = "linreg", "formula" = stanread7 ~ read2 + hispanic + read9 + read1.cgm + lrnprob1.cgm + behsymp1, nodes = nodes.stanread7)
# combine predictor models into a list
predictor.models <- list(read2 = model.read2, hispanic = model.hispanic, read9 = model.read9, behsymp1 = model.behsymp1, lrnprob1 = model.lrnprob1, read1 = model.read1)
# estimate factored regression model w mdmb
fit <- frm_em(dat = data, dep = model.stanread7, ind = predictor.models)
summary(fit)
source("~/Library/CloudStorage/Dropbox/Project Folder (1)/Analysis Scripts/Ex4. FIML Moderated Regression/Ex4.R")
source("~/Library/CloudStorage/Dropbox/Project Folder (1)/Analysis Scripts/Ex4. FIML Moderated Regression/Ex4.R")
# Three-Level Regression with a Cross-Level Interaction
library(fdir)
library(mitml)
library(lme4)
# FIML Moderated Regression
library(fdir)
library(lavaan)
library(mdmb)
##################################################
# read data
##################################################
# set working directory
set()
# read raw data from working directory
data <- read.table("behaviorachievement.dat", na.strings = "999")
names(data) <- c("id","male","hispanic","riskgrp","atrisk","behsymp1","lrnprob1",
"read1","read2","read3","read9","readgrp9","stanread7",
"math1","math2","math3","math9","mathgrp9","stanmath7")
##################################################
# estimate means and center variables
##################################################
# estimate sample statistics in lavaan
model <- "stanread7 ~ 1; read2 ~ 1; hispanic ~ 1; read9 ~ 1; read1 ~ 1; lrnprob1 ~ 1; atrisk ~ 1;"
descriptives <- inspectSampleCov(model, data, missing = "fiml")
# center lower-order variables
data$read1.cgm <- data$read1 - descriptives$mean["read1"]
data$lrnprob1.cgm <- data$lrnprob1 - descriptives$mean["lrnprob1"]
##################################################
# analysis
##################################################
# summarize incomplete predictors to determine ranges for pseudo-imputations
summary(data[,c("stanread7","read2","read9","read1.cgm","lrnprob1.cgm","atrisk")])
# set ranges (nodes) for pseudo-imputations
nodes.stanread7 <- seq(80, 420, by = 5)
nodes.read2 <- seq(10, 160, by = 2)
nodes.read9 <- seq(30, 130, by = 2)
nodes.read1 <- seq(-55, 75, by = 2)
nodes.lrnprob1 <- seq(-30, 50, by = 2)
nodes.atrisk <- c(0,1)
# model for behsymp1 predictor
model.atrisk <- list( "model" = "logistic", "formula" = atrisk ~ 1, nodes = nodes.atrisk)
# model for lrnprob1 predictor
model.lrnprob1 <- list( "model" = "linreg", "formula" = lrnprob1.cgm ~ atrisk, nodes = nodes.lrnprob1)
# model for read1 predictor
model.read1 <- list( "model" = "linreg", "formula" = read1.cgm ~ lrnprob1.cgm + atrisk, nodes = nodes.read1)
# focal model for read9 outcome
model.read9 <- list( "model" = "linreg", "formula" = read9 ~ read1.cgm + lrnprob1.cgm + read1.cgm*lrnprob1.cgm + atrisk, nodes = nodes.read9)
# model for read2 auxiliary variable
model.read2 <- list("model" = "linreg", "formula" = read2 ~ read9 + read1.cgm + lrnprob1.cgm + atrisk, nodes = nodes.read2)
# combine all previous models into a list
predictor.models <- list(atrisk = model.atrisk, lrnprob1 = model.lrnprob1, read1 = model.read1, read9 = model.read9, read2 = model.read2)
# final model for stanread7 auxiliary variable
model.stanread7 <- list("model" = "linreg", "formula" = stanread7 ~ read2 + read9 + read1.cgm + lrnprob1.cgm + atrisk, nodes = nodes.stanread7)
# estimate factored regression model w mdmb
fit <- frm_em(dat = data, dep = model.stanread7, ind = predictor.models)
summary(fit)
# FIML Moderated Regression
library(fdir)
library(lavaan)
library(mdmb)
##################################################
# read data
##################################################
# set working directory
set()
# read raw data from working directory
data <- read.table("behaviorachievement.dat", na.strings = "999")
names(data) <- c("id","male","hispanic","riskgrp","atrisk","behsymp1","lrnprob1",
"read1","read2","read3","read9","readgrp9","stanread7",
"math1","math2","math3","math9","mathgrp9","stanmath7")
##################################################
# estimate means and center variables
##################################################
# estimate sample statistics in lavaan
model <- "stanread7 ~ 1; read2 ~ 1; read9 ~ 1; read1 ~ 1; lrnprob1 ~ 1; atrisk ~ 1;"
descriptives <- inspectSampleCov(model, data, missing = "fiml")
# center lower-order variables
data$read1.cgm <- data$read1 - descriptives$mean["read1"]
data$lrnprob1.cgm <- data$lrnprob1 - descriptives$mean["lrnprob1"]
##################################################
# analysis
##################################################
# summarize incomplete predictors to determine ranges for pseudo-imputations
summary(data[,c("stanread7","read2","read9","read1.cgm","lrnprob1.cgm","atrisk")])
# set ranges (nodes) for pseudo-imputations
nodes.stanread7 <- seq(80, 420, by = 5)
nodes.read2 <- seq(10, 160, by = 2)
nodes.read9 <- seq(30, 130, by = 2)
nodes.read1 <- seq(-55, 75, by = 2)
nodes.lrnprob1 <- seq(-30, 50, by = 2)
nodes.atrisk <- c(0,1)
# model for behsymp1 predictor
model.atrisk <- list( "model" = "logistic", "formula" = atrisk ~ 1, nodes = nodes.atrisk)
# model for lrnprob1 predictor
model.lrnprob1 <- list( "model" = "linreg", "formula" = lrnprob1.cgm ~ atrisk, nodes = nodes.lrnprob1)
# model for read1 predictor
model.read1 <- list( "model" = "linreg", "formula" = read1.cgm ~ lrnprob1.cgm + atrisk, nodes = nodes.read1)
# focal model for read9 outcome
model.read9 <- list( "model" = "linreg", "formula" = read9 ~ read1.cgm + lrnprob1.cgm + read1.cgm*lrnprob1.cgm + atrisk, nodes = nodes.read9)
# model for read2 auxiliary variable
model.read2 <- list("model" = "linreg", "formula" = read2 ~ read9 + read1.cgm + lrnprob1.cgm + atrisk, nodes = nodes.read2)
# model for stanread7 auxiliary variable
model.stanread7 <- list("model" = "linreg", "formula" = stanread7 ~ read2 + read9 + read1.cgm + lrnprob1.cgm + atrisk, nodes = nodes.stanread7)
# combine predictor models into a list
predictor.models <- list(atrisk = model.atrisk, lrnprob1 = model.lrnprob1, read1 = model.read1, read9 = model.read9, read2 = model.read2)
# estimate factored regression model w mdmb
fit <- frm_em(dat = data, dep = model.stanread7, ind = predictor.models)
summary(fit)
# summarize incomplete predictors to determine ranges for pseudo-imputations
summary(data[,c("stanread7","read2","read9","read1.cgm","lrnprob1.cgm","atrisk")])
